function [] = calcB0(subjectID, b0ID)
%calcB0 Calculate B0 from gre data with the phase-difference method using two TE data points.
%  Filenames are are subjectID_b0ID_mri.mnc or subjectID_b0ID(1)_mri.mnc
%  and subjectID_b0ID(2)_mri.mnc
%
%  Author: Mathieu Boudreau
%  Date: August 21 2013

  
    % Split B0 phase images
    if(length(b0ID)==1)
        system(sprintf('mincreshape -clobber  -dimrange time=0,0 %s_%i_mri.mnc  b0_te_1.mnc', subjectID, b0ID))
        system(sprintf('mincreshape -clobber  -dimrange time=1,0 %s_%i_mri.mnc  b0_te_2.mnc', subjectID, b0ID))
    % Copy and rename both B0 phase files if two files were inputed.
    elseif((length(b0ID)==2))
        system(sprintf('cp %s_%s_mri.mnc b0_te_1.mnc',subjectID, num2str(b0ID(1))))
        system(sprintf('cp %s_%s_mri.mnc b0_te_2.mnc',subjectID, num2str(b0ID(2))))
    else
        disp('qmt_preprocess can only process gre_b0 B0 data. Stop.')
        return
    end
    
    calculateB0FieldMap_3t('b0_te_1.mnc', 'b0_te_2.mnc', 'b0_map.mnc')
    system('mv b0_map.mnc b0/b0_map.mnc')

end

